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responses of tropical plant to cold stress 
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Abstract 

Background: Cassava is an important tropical root crop adapted to a wide range of environmental stimuli such as 
drought and acid soils. Nevertheless, it is an extremely cold-sensitive tropical species. Thus far, there is limited 
information about gene regulation and signalling pathways related to the cold stress response in cassava. The 
development of microarray technology has accelerated the study of global transcription profiling under certain 
conditions. 

Results: A 60-mer oligonucleotide microarray representing 20,840 genes was used to perform transcriptome 
profiling in apical shoots of cassava subjected to cold at 7°C for 0, 4 and 9 h. A total of 508 transcripts were 
identified as early cold-responsive genes in which 319 sequences had functional descriptions when aligned with 
Arobidopsis proteins. Gene ontology annotation analysis identified many cold-relevant categories, including 
'Response to abiotic and biotic stimulus', 'Response to stress', 'Transcription factor activity', and 'Chloroplast'. Various 
stress-associated genes with a wide range of biological functions were found, such as signal transduction 
components (e.g., MAP kinase 4), transcription factors (TFs, e.g., RAP2.n), and reactive oxygen species (ROS) 
scavenging enzymes (e.g., catalase 2), as well as photosynthesis-related genes (e.g., PsoL). Seventeen major TF 
families including many well-studied members (e.g., AP2-EREBP) were also involved in the early response to cold 
stress. Meanwhile, KEGG pathway analysis uncovered many important pathways, such as 'Plant hormone signal 
transduction' and 'Starch and sucrose metabolism'. Furthermore, the expression changes of 32 genes under cold 
and other abiotic stress conditions were validated by real-time RT-PCR. Importantly, most of the tested stress- 
responsive genes were primarily expressed in mature leaves, stem cambia, and fibrous roots rather than apical 
buds and young leaves. As a response to cold stress in cassava, an increase in transcripts and enzyme activities of 
ROS scavenging genes and the accumulation of total soluble sugars (including sucrose and glucose) were also 
detected. 

Conclusions: The dynamic expression changes reflect the integrative controlling and transcriptome regulation of 
the networks in the cold stress response of cassava. The biological processes involved in the signal perception and 
physiological response might shed light on the molecular mechanisms related to cold tolerance in tropical plants 
and provide useful candidate genes for genetic improvement. 
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Background 

Cassava {Manihot esculenta Crantz) is widely cultivated 
for its starchy storage roots and is a staple food and ani- 
mal feed in tropical and sub-tropical areas [1]. It is also 
considered to be an important source of modified 
starches and bioethanol in China and other Southeast 
Asian countries [2,3]. Nevertheless, as a tropical root 
crop, cassava is native to a warm habitat and is categor- 
ized as a cold-sensitive species [4]. Thus, low tempera- 
tures and frozen conditions are the most important 
limiting factors for its geographical location and produc- 
tivity. In the subtropics, where unpredictable cold 
weather occurs occasionally, it is important to protect 
the storage roots and propagation stems from chilling 
stress. For example, the unprecedented freezing disaster 
occurred in Southern China in January 2008 caused 
great damage to cassava stem seeds and led to yield 
reduction in Guangxi, Guangdong and other provinces, 
resulting in a loss of a billion Chinese Yuan [5]. Further- 
more, to ensure a prolonged growth period (i.e., early 
planting and late harvesting) in the high latitude regions, 
novel cassava cultivars with improved cold tolerance are 
in demand. 

Under low temperature below 10°C, many species of 
tropical or subtropical origin are typically injured or 
killed and show various symptoms of chilling injury due 
to the inability to adapt to non-freezing low tempera- 
tures [6]. For example, cassava exhibits obvious symp- 
toms of damage at these temperatures, including 
delayed sprouting of the stem cutting, yield decrease, 
reduced leaf expansion, chlorosis and even necrosis in 
its leaves [7]. Low temperatures have also been recog- 
nized as an important facilitator of decreases in nutrient 
absorption rates (e.g.. Boron), reductions in the leaf 
photosynthetic rate, and the inhibition of plant growth 
[4]. In addition, the physiological status of cold-stressed 
plants is also altered, such as transient increases in hor- 
mone levels (e.g., ABA) [8] and changes in membrane 
lipid composition [9]. Furthermore, the accumulation of 
compatible osmolytes, such as soluble sugars, betaine, 
and proline [10-12], and increases in the level of antioxi- 
dants [13] are also occurred. In contrast, temperate 
plants can withstand freezing temperatures following a 
period of low, but non-freezing, temperatures, a process 
called cold acclimation [14]. The mechanisms of cold 
acclimation have been extensively investigated in the 
model plant Arabidopsis thaliana [15], as well as in 
other important crop species such as maize and barley 
[16,17]. The most extensively studied mechanism 
involves a class of ethylene response factors (ERFs) 
known as the dehydration-responsive element-binding 
proteins/C-repeat-binding factors (DREBs/CBFs). They 
interact with the dehydration-responsive element/C- 
repeat element (DRE/CRT) in the promoters of their 



downstream target genes to execute a highly coordi- 
nated transcriptional response to low temperature sig- 
nals. Overall, less information is available for tropical 
plants (e.g., cassava). Therefore, it would be important 
to increase our understanding of the physiological, bio- 
chemical, and molecular characteristics that are asso- 
ciated with cold stress in tropical cassava, especially 
during the early stages. Moreover, the function of 
homologous genes in Arabidopsis could be used to pre- 
dict the function of cassava genes which may share the 
same biological functions or common regulatory scenar- 
ios among different plant species [18]. 

With the development of molecular technologies and 
'omics' tools, microarray studies have become a useful 
strategy for the global analysis of plant gene expression. 
Using cDNA microarrays or whole genome arrays, the 
abiotic stress responses of Arabidopsis and other plants 
have been widely analyzed. For examples, the expression 
patterns of genes were identified in Arabidopsis and rice 
under conditions of drought, cold, high-salinity or absci- 
sic acid treatment [19,20]. The transcriptomic identifica- 
tion of candidate genes involved in sunflower responses 
to chilling and salt stresses has also been reported [21]. 
The long-oligonucleotide microarray, which has many 
advantages compared to other microarray technologies, 
reliably detects transcript ratios even at one copy per 
cell in complicated biological samples, and could distin- 
guish different members of gene families [22,23]. Thus, 
it is appropriate to explore the early cold stress respon- 
sive genes of cassava using a 60-mer long-oligonucleo- 
tide microarray approach. 

Thus far, several studies have been reported in cassava 
using genomic tools for developing new knowledge and 
technologies. For instance, EST and cDNA libraries have 
been constructed in cassava for the identification of 
starch biosynthesis and biotic/abiotic-responsive genes 
or genes corresponding to the different developmental 
stages of various tissues [24-28]. Studies using cassava 
cDNA microarrays to analyze gene expression in cassava 
plants subjected to post-harvest physiological deteriora- 
tion (PPD) or Xanthomonas axonopodis infection have 
also been reported [29,30]. Moreover, gene expression 
profiling related to the different growth stages of cassava 
storage roots was recently conducted using a long-oligo- 
nucleotide microarray [31]. Fortunately, a draft genome 
sequence for cassava has been released and recently 
updated (http://www.phytozome.net/cassava), which 
greatly facilitates cassava research worldwide. Neverthe- 
less, comprehensive genome-wide expression profiling 
data for cassava under cold stress treatment are still 
lacking, which limits our capacity to decipher the mole- 
cular mechanisms related to various stress responses. 
Additionally, to stabilize cassava yield even under unfa- 
vorable growth conditions, it is necessary to study the 
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changes of pan-genome genes or pathways in order to 
identify key candidates for improving cold tolerance in 
cassava. 

In this study, a gene expression profiling of apical 
shoots from cassava subjected to low temperature was 
conducted using a custom-designed, 60-mer oligonu- 
cleotide microarray covering 20,840 cassava transcripts. 
In total, 508 differentially expressed genes (DEGs) were 
identified, and their biological functions were character- 
ized through gene ontology (GO) annotation and KEGG 
pathway analysis. Our study could increase the under- 
standing of cassava gene regulation under the condition 
of cold stress and reveal novel approaches to improve 
cold tolerance through genetic engineering. 

Results and discussion 

Phenotypic and physiological changes of cassava in 
response to cold 

Cassava is a typical tropical crop adapted to warm cli- 
mates. Under normal conditions, 3-month-old cassava 
(TMS60444) plants have vigorous apical buds and fully 
expanded mature leaves with upward leaf fingers (Figure 
lA, Left panel). When a low-temperature treatment is 
applied to cassava at 7°C under weak light for 4 h, the 
plants displayed visible morphological changes, includ- 
ing weak dehydration and wilting of the apical buds and 
leaves, especially in the buds, as well as downward leaf 
fingers (Figure lA, Middle panel). Under prolonged 
exposure to the stress (9 h), the whole plants exhibited 
obvious phenotypic damages, including softening and 
downward bending of the petioles, loss of strength in 
the immature stems and more severe wilting leaves (Fig- 
ure lA, Right panel). Generally, the treated plants could 
be partially recovered and resumed growth when trans- 
ferred into normal conditions (data not shown), indicat- 
ing that the cellular changes caused by a 9 h cold 
treatment were reversible. Since apical shoots show 
extremely sensitive response to low temperature (Addi- 
tional file 1), the apical shoots (about 4-6 cm) of 3- 
month-old plants, which were treated for 0, 4 and 9 h, 
were used as material for microarray hybridization. 

In addition to morphological changes, a clear impact 
on ultra-structure of chloroplasts was also observed in 
leaves of cold-stressed cassava plants. Compared to nor- 
mal chloroplasts in unstressed leaves (Figure IB), the 
stacking numbers of thylakoids were dramatically 
reduced and less organized, and the starch granules had 
almost disappeared (Figure IC). The ultra-structural 
changes in cellular organelles indicate that cold stress 
might exert an adverse impact on plant growth, causing 
a reduced photosynthetic rate under the stress. 

Cold stress also leads to obviously physiological 
changes in cassava. Two stress responsive metabolites, 
namely malondialdehyde (MDA) and proline, were 



monitored. MDA is considered to be the final product 
of lipid peroxidation in the plant cell membrane and is 
an important indicator of membrane system injuries and 
cellular metabolism deterioration [32]. It has also been 
reported that MDA contents change during cold stress 
treatment in plants [33]. In our experiment, compared 
to the control level, the MDA concentrations decreased 
rapidly to 50% after 4 h but increased almost 25% after 
9 h; the observed levels had experienced a nearly two- 
fold change after 24 h of stress treatment (Figure ID). 
Decrease in MDA content at 4 h might be due to a 
transient stress response of cassava to the low tempera- 
ture shock. But, the prolonged treatment (9 h) finally 
led to cell damages and MDA accumulation. The accu- 
mulation of proline is frequently associated with whole 
plant tolerance to chilling and other stresses [34]. Unlike 
that of MDA, the concentration of proline increased 
rapidly and achieved nearly a three-fold change after 4 
h; it then decreased slightly but remained at a high level 
until 24 h (Figure IE). These results were consistent 
with previous studies that proline accumulated in leaves 
exposure to cold, salt, and other stresses in Arabidopsis 
and other plants [35]. The accumulation of these meta- 
bolites is a good indication that the stress-treated plants 
were actively mounting a stress response during the per- 
iods when they were subjected to cold stress. 

Microarray for transcriptomic analysis of low-temperature 
treated cassava 

To explore the transcriptomic changes of cassava in 
response to cold stress, a custom long-oligonucleotide 
(60-mer) microarray generated by the Agilent SurePrint 
ink-jet technology was applied in the current study. The 
detailed protocol for designing the cassava microarray has 
been described by Yang et al [31] To validate the microar- 
ray quality, the signal-to-noise ratio (SNR) for each spot 
was calculated as described by Leiske et al [23] The aver- 
age percentages of acceptable spots (SNR > 2.6) and high- 
quality spots (SNR > 10) were 94.53 ± 1.84% and 79.24 ± 
1.76% in all 9 arrays, which demonstrated the overall 
reproducibility and high quality of the array. The hierarch- 
ical clustering results indicated the clear separation of the 
samples from different time points (Figure 2A), suggesting 
that the entire experiment from sample collection to data 
extraction was reproducible and reliable. Five house-keep- 
ing genes (beta-actin, cl5, EFla, RuBisCO small chain pre- 
cursor, and TUB6) were used as the internal controls and 
their expressions were consistent in the samples at differ- 
ent time points (Additional file 2), confirming the reliabil- 
ity and accuracy of these microarrays. 

Transcriptomic responses to cold treatment 

Before normalization, the signal intensities of each fea- 
ture were filtered against negative controls on the array. 
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Figure 1 Phenotypic and physiological changes in cold-stressed cassava. (A) 3-month-old cassava subjected to low temperature stress (7°C) 
for 0, 4, and 9 h in a chamber under weak light showing phenotypic changes. (B) Fully expanded cassava leaf (Left panel) and TEM analysis of a 
chloroplast (CP1, Right panel). (C) Low-temperature treated cassava leaf showing dehydration and necrosis (Left panel) and TEM analysis showing 
chloroplast abnormalities (CP2, Right panel). (D) MDA contents in cassava apical leaves exposed to 7°C for 0, 4, 9, and 24 h. (E) Proline 
accumulation in cassava apical leaves exposed to 7°C for 0, 4, 9, and 24 h. The double asterisks indicate a statistically significant difference (p < 
0.01) for the data of the stress-treated samples compared to those of the unstressed samples. The mean values are calculated from three 
biological replicates; the error bars represent the standard error of the mean (SEM). Bar = 500 nm. 
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Figure 2 Expression profiling of cold-regulated genes in cassava apical shoots. (A) Hierarchical cluster analysis. (B) Venn diagrams showing 
cold-regulated genes across three comparisons (4 h/0 h, 9 h/0 h, and 9 h/4 h). The red numbers are the total numbers of differentially 
expressed genes (DEGs); the percentages in parentheses were calculated as the ratio of regulated genes to the total number of cold-regulated 
genes (508). It should be noted that one gene was up-regulated at 4 h/0 h but down-regulated at 9 h/4 h on our array. 
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and 69.87 ± 0.72%, 68.14 ± 2.59%, and 71.50 ± 5.75% of 
probes on the array were expressed at 0, 4 and 9 h, 
respectively. Two filtering criteria were used to define 
differentially expressed genes (DEGs) in our data analy- 
sis: a two-fold change in transcript levels among every 
two time points and a P value < 0.05. To analyze the 
similarities and differences among the cold-responsive 
transcriptomes, a hierarchical clustering was prepared to 
represent the transcripts of all the differentially 
expressed probes at the 3 replicates of 0, 4 and 9 h. 
These results indicated significant differences in the 
gene expression profiles between the treatments of 9 h 
and 0 h or 9 h and 4 h, in contrast to the relatively high 
similarity of the expression profiles between 0 h and 4 h 
(Figure 2A). Among the DEGs, 58, 252, and 92 genes 
were up-regulated and 13, 210, and 36 were down-regu- 
lated by analyzing 4 h/0 h, 9 h/0 h, and 9 h/4 h, respec- 
tively (Figure 2B). Notably, compared to only 71 and 
128 genes showed differentially expressed at 4 h/0 h and 
9 h/4 h, 474 DEGs were found at 9 h/0 h, suggesting a 
prolonged stress treatment (9 h) could trigger more 
stress-related gene expression. Moreover, about 80% of 
4 h/0 h and 75% of 9 h/4 h DEGs overlapped with those 
of 9 h/0 h, which indicated a strong linkage among the 
3 stressed comparison points and a progressive biologi- 
cal process. 

In total, 508 DEGs, including 292 cold-inducible and 
215 cold-repressed genes, were identified on our array; 
only one gene (CUST_14887) was up-regulated at 4 h/0 
h but down-regulated at 9 h/4 h (Additional file 3). 
These represent about 3.5% of the total expressed cas- 
sava transcripts on the array. This result indicated that 
even for tropical crops (e.g., cassava), which have 
adapted to warm climates, possess a large amount of 
cold-responsive genes as other temperate plants (e.g., 
Arabidopsis), Therefore, the differences in the ability of 
tropical species and temperate species to tolerate cold 
might not totally be due to the amount of genes respon- 
sive to cold stress. Instead, we speculate that plants 
require to timely and coordinately mobilize the biologi- 
cal functions and regulatory networks of stress-respon- 
sive genes against stress. Consistent with our 
observations, a recent study confirmed such scenario 
even between two close plant species [36]. 

Gene ontology clustering of cold-regulated genes 

To explore the biological functions of cold-responsive 
genes, a total of 20,840 sequences were used as a query 
to perform an alignment with Arabidopsis proteins. Of 
these, 20,450 sequences had hits with 11,995 Arabidop- 
sis proteins, and 14,813 (about 71.8%) sequences had 
hits with 9,158 Arabidopsis proteins at an E value < le- 
5. After searching the descriptions of the 508 DEGs, 319 
queries had the highest homologies with annotated 



proteins according to The Arabidopsis Information 
Resource (TAIR), which accounted for 62.8% of the 
responsive genes (Additional file 4). To determine the 
detailed function of the 319 hits with Arabidopsis gene 
locus identifiers, a GO annotation was performed with 
the GO terms of TAIR (http://www.Arabidopsis.org/ 
tools/bulk/go/index.jsp), in which 370 GO IDs were hit 
by 291 gene locus identifiers (Additional file 5). There- 
fore, functional clusters that were classified according to 
the three components ('Biological Process', 'Molecular 
Function', and 'Cellular Component') are presented and 
discussed (Additional files 6, 7, and 8, respectively). 
Based on the TAIR percent analysis, the categories of 
'Response to abiotic and biotic stimulus' and 'Response 
to stress', 'Transcription factor activity', and 'Chloro- 
plast' were the major cold-related GO Slims that 
strongly attracted our attention. 

Cold-responsive genes related to 'Response to abiotic and 
biotic stimulus' and 'Response to stress' 

In the GO terms of 'Biological Process', the categories of 
'Response to abiotic and biotic stimulus' and 'Response 
to stress', both accounted for 13.36% of the total GOs, 
were suggested to be the most relevant to cold stress 
(Figure 3A). The largest proportions of 'Biological Pro- 
cess' terms were annotated as 'Other cellular processes' 
(48.09%), 'Other metabolic processes' (43.51%), 
'Unknown biological process' (33.59%), and 'Other biolo- 
gical processes' (14.12%), indicating comprehensive 
changes in cassava gene expression (Figure 3A). In addi- 
tion, 'Protein metabolism' was another noticeable major 
category (15.27%) supposed to be also involved in the 
cold stress response [37]. 

There were 34 and 33 genes assigned to the GO terms 
of 'Response to abiotic and biotic stimulus' and 
'Response to stress', respectively (Additional file 9). 
Because 27 genes were common in both categories, we 
combined the two categories for analysis and discussion 
(Table 1). Among the overlapping genes, 20 up-regu- 
lated genes and 7 down-regulated genes were included. 
Their encoding proteins cover a wide range of biological 
functions, including transcription factor (TFs, e.g., AP2 
domain transcription factor), protein kinases (e.g., pro- 
Une extension-like receptor kinase 1), small heat shock 
proteins (e.g., heat shock protein 18.2), ROS scavenging 
enzymes [e.g., catalase 2 {CAT2)], fatty acid desaturases 
[e.g., fatty acid biosynthesis 2 {SSI2)], signaling molecu- 
lar proteins (e.g., MAP kinase 4), programmed cell 
death (PCD)-related gene (e.g., chitinase), and genes 
involved in the gibberelUn acid (GA) and jasmonic acid 
(JA) metabolism pathway [e.g., gibberellic acid insensi- 
tive (GA/)]. Among the genes specifically related to the 
'Response to abiotic and biotic stimulus', protein degra- 
dation genes (e.g., ubiquitin-protein ligase) and early 
flowering 4 were up-regulated. Meanwhile, genes 
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Figure 3 TAIR percent of gene ontology (GO) terms for (A) 'Biological Process', (B) 'Molecular Function', and (C) 'Cellular Component' 
of cold up-regulated and down-regulated genes. 



involved in ABA signal transduction [e.g., ABA-respon- 
sive element binding protein 3 (AREB3)] and disease 
resistance protein (CC-NBS-LRR), which were down- 
regulated by cold, were uniquely identified in GO Slim 



as 'Response to stress' (Table 1). Most of these cold- 
responsive genes were either up-regulated or down- 
regulated at 9 h/0 h, except 4 genes encoding the cal- 
cium ion binding protein (TCH2), ethylene response 
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Table 1 Low temperature responsive transcripts related to 'Response to abiotic and biotic stimulus' and 'Response to 
stress' 



Probe Name 


Description 


4 h/0 h 


9 h/0 h 
Log2 Ratio 


9 h/4 h 


AG! Locus 


E Value 


'Response to abiotic and biotic stimulus' and 'Response to stress' 












r"l ICT IQQO 


AP2 domain transcription factor 




4.02 




AT4G34410 


1 .60E-23 


r-\ 1 CT 1 /I 1 1 o 

LUb l_ 14 1 1 o 


SSI2 (fatty acid biosynthesis 2) 


2.31 


3.69 


1.38 


AT2G43710 


1.80E-57 


r-\ ICT 1 


alpha-crystallin domain 31.2 




2.73 


2.33 


AT1G06460 


5.80E-21 


r-\ ic~r 1 moo 


MAP l<inase 4 




2.55 




AT4G01370 


1.30E-16 


r-\ ic~r cicn 


proline extensin-lil<e receptor kinase 1 


1.18 


2.46 


1.28 


AT3G24550 


3.40E-29 


r-\ ic~r 1 1 1 

LUb l_l 1 1 16 


heat shocl< protein 18.2 




2.23 


1.59 


AT5G59720 


2.10E-48 


LUb l_oo 1 U 


Bet V 1 allergen family protein 


1.81 


2.16 




AT5G28010 


4.90E-06 


CUST_7667 


chitinase 




2.04 




AT3G 12500 


1.50E-11 


^\ 1 c — r 1 1 

CUbl_l 1923 


MEE14 (maternal effect embryo arrest 14) 




1.85 


1.81 


AT2G15890 


3.40E-42 


r~\ ICT ^C^O 

LUb l_Dbob 


MLP-like protein 31 




1.72 




AT1G70840 


1 .OOE-43 


/^i 1 CT 1 cn/i ^ 

LUb l_ 1 by4D 


MLP-like protein 31 




1.42 




AT1G70840 


2.30E-38 


r~| ICT 1 1 QHQ 

LUbl_l IbUb 


MLP-like protein 34 




1.42 




AT1G70850 


2.60E-24 


r"l ICT 0 CO /! 

LUb I_bbb4 


Jasmonic acid carboxyl methyltransferase 




1.18 




ATI G 19640 


5.40E-40 


r"! 1 CT 1 r\nno 

LUb l_ 1 Uyyo 


MLP-like protein 423 




1.14 




AT1G24020 


2.00E-10 


r'\ 1 CT 1 1^ CO 
LUb l_ 1 zoby 


beta-ketoacyl-CoA synthase 




1.11 




AT5G43760 


1 .40E-63 


CUST_3051 


CAT2 (catalase 2) 




1.11 


1.01 


AT4G35090 


7.60E-76 


CUST_14725 


heat shock protein 17.4 




1.06 




AT3G46230 


4.10E-33 


CUST_17985 


TCH2 (calcium ion binding protein) 






1.09 


AT5G37770 


2.10E-39 


CUST_14955 


ERF2 (ethylene response transcription activator 2) 






3.35 


AT5G47220 


1.50E-33 


CUST_7141 


MAP kinase kinase 9 






1.46 


AT1G73500 


1.70E-42 


CUST_14327 


RING-H2 finger A2A zinc ion binding 




-1.04 




ATI Gl 51 00 


5.40E-14 


CUST_13235 


GAI (gibberellic acid insensitive) transcription factor 




-1.07 




ATI G 14920 


2.00E-108 


CUST_14166 


GSTU19 (glutathione transferase) 




-1.12 




AT1G78380 


1.90E-70 


CUST_9208 


UDP-glucosyltransferase 




-1.25 




AT4G01070 


8.80E-42 


CUST_6759 


MYB family transcription factor 


_ 


-1.49 


_ 


AT1G49010 


9.30E-41 


CUST_5609 


RING-H2 finger A2A zinc ion binding 




-1.97 




ATI Gl 51 00 


6.70E-39 


CUST_18370 


CBL-interacting protein kinase 25 




-2.03 




AT5G25110 


9.90E-17 


'Response to abiotic and biotic stimulus' 












CUST_2683 


early flowering 4 


_ 


4.01 


3.64 


AT2G40080 


1 .OOE-25 


CUST_4383 


gibberellin 2-beta-dioxygenase 




3.63 




AT1G30040 


4.70E-51 


CUST_11786 


ubiquitin-protein ligase 




3.48 


3.32 


AT1G68050 


8.90E-60 


CUST_12768 


pesudo-response regulator 5 






1.31 


AT5G24470 


2.80E-20 


CUST_14809 


SIGB (SIGMA FACTOR B) 


- 


-1.03 




AT1G08540 


5.30E-20 


CUST_14095 


polyubiquitin 3 protein binding 




-1.05 




AT5G03240 


8.90E-17 


CUST_11479 


zinc ion binding 




-1.18 


-1.01 


AT1G52520 


1 .OOE-72 


'Response to stress' 












CUST_8448 


peroxidase 72 




3.42 


2.57 


AT5G66390 


3.90E-35 


CUST_8771 


Homolog of yeast autophagy 18 




2.24 




AT3G62770 


2.80E-63 


CUST_11592 


ozone-responsive stress-related protein 




1.72 




AT1G01170 


7.10E-27 


CUST_10706 


ABA-responsive element binding protein 3 




-1.18 




AT3G56850 


1 .20E-35 


CUST_11310 


disease resistance protein (CC-NBS-LRR) 




-1.18 




AT1G61190 


6.40E-06 


CUST_18065 


ABA-responsive element binding protein 3 




-1.22 




AT3G56850 


2.30E-20 



Notes: Low temperature responsive genes are grouped to 'Response to abiotic and biotic stimulus' and 'Response to stress' 
and negative values of Log2 Ratio are either up- or down-regulated genes In the three pair-comparisons. No significant fold 



based on the GO annotation. Positive 
changes are Indicated by 



transcription activator 2 (ERF2), MAP kinase kinase 9, 
and pseudo-response regulator 5. It is notable that 27 
out of these 40 stress-associated genes were up-regu- 
lated, indicating that the up-regulation of cold- 



responsive transcripts by low temperature may play a 
crucial role in the stress tolerance of plants [38]. 
Of their homologous genes in Arabidopsis, the MAP 
kinase 4, beta-ketoacyl-CoA synthase, CAT2, and TCH2 
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were previously reported to be involved in the cold 
stress response process. For example, MAP kinase 4, 
which is required for cytokinesis, was implicated in cold 
and salt stress tolerance [39]; beta-ketoacyl-CoA 
synthase, which is involved in the biosynthesis of 
VLCFA (very long chain fatty acids), was also responsive 
to cold and other osmotic stresses [40]. CAT2, a target 
of plastid, could increase its transcript abundance in 
response to cold stress [41]. TCH2 has been implicated 
in calcium signaling in response to diverse stimuli, such 
as cold, light, pathogens, and touch [42]. Gibberellin 2- 
beta-dioxygenase, which degrades active GAs, is 
involved in cold response [43] and up-regulation of its 
transcript implied that GA level may be reduced under 
cold stress. Additionally, calcineurin B-like (CBL) -inter- 
acting protein kinase, which has protein serine/threo- 
nine kinase activity, has been reported to be involved in 
the defense response to fungus [44]. These results 
revealed that there are considerable conserved and com- 
mon components in cold stress response mechanisms 
across plant species, including temperate and tropical 
plants. 

Besides, a large number of genes with known function 
were firstly identified in the cold response. The AP2 
domain transcription factor (CUST_2332), a member of 
the ERF (ethylene response factor) subfamily B-3 of 
ERF/AP2 transcription factor family, has been reported 
involved in response to fungus and chitin in Arabidop- 
sis. In this study, the transcript abundance of the AP2 
domain transcription factor was up-regulated 16-fold by 
cold, indicating that it may play a significant role in the 
regulation of novel signaling pathways to enhance cas- 
sava cold tolerance. The expression of a stearoyl-ACP 
desaturase gene SSI2, which involves in fatty acid desa- 
turation, was up-regulated almost 10-fold under cold 
stress. This change suggested that the gene might be 
participated in altering membrane lipid composition to 
enhance membrane fluidity. Maternal effect embryo 
arrest 14 {MEE14), which required for embryonic devel- 
opment ending in seed dormancy, was also highly 
induced in this study. Interestingly, three transcripts 
(early flowering 4, gibberellin 2-beta-dioxygenase, and 
pseudo-response regulator 5) that respond to red or far- 
red light in Arabidopsis were cold-inducible, indicating a 
crosstalk between light signaling and cold stress 
response. 

In addition, a large number of genes with unknown 
function were also firstly recorded as cold responsive 
genes. For instance, four members of MLP-like protein 
(CUST_6563, CUST_15946, CUST_11303, and 
CUST_10998) were all up-regulated by 9 h cold treat- 
ment. Alpha-crystallin domain 31.2, Bet v I allergen 
family protein, and homology of yeast autophagy 18 
were highly induced more than 5-fold by cold. Two 



members of RING-H2 finger A2A zinc ion binding pro- 
teins were also found to be down-regulated. Taken 
together, cassava might evolve not only conserved but 
also specific molecular mechanisms related to stress sig- 
naling and response. 

Transcription factors response to cold stress 

For the 'Molecular Function' GO terms. Transferase 
activity', 'Hydrolase activity', and 'Transcription factor 
activity' were the three major categories. They 
accounted for 13.28%, 11.72%, and 11.33% of the total 
GOs, respectively (Figure 3B). Among them, genes asso- 
ciated with 'Transcription factor activity' may be central 
regulators involved in early cold signal transduction that 
trigger a cascade of downstream gene expression. 

Transcription factors (TFs) play a significant role in 
plant development and stress tolerance [38]. To identify 
the transcription factors involved in the cold stress 
response, we surveyed the biological functions of puta- 
tive TFs that were differentially expressed in cassava 
under the cold treatment. A total of 32 genes were iden- 
tified as TFs when the 319 identifiers were annotated in 
TAIR. The number of up-regulated TFs (15) was nearly 
equal to the number of down-regulated ones (17), sug- 
gesting that both transcriptional activation and repres- 
sion are involved. The Arabidopsis genome has more 
than 2,657 predicted TFs belonging to 81 gene families 
in the Plant Transcription Factor Database (PlnTFDB, 
http://plntfdb.bio. uni-potsdam.de/v3. 0/index.php?sp_i- 
d=ATH) [45]. According to this classification, 31 cold- 
responsive transcription factors fell into 17 families, 
except 1 (CUST_19860) with non-classification (Table 
2), and they were annotated with the detailed GO Slims 
(Additional file 10). In Arabidopsis, at least five TF 
families have been reported to be involved in the cold 
stress response process, including AP2-EREBP (APE- 
TALA2/ET-Responsive Element Binding Protein), MYB 
(Myeloblastosis), NAC (NAM, ATAFl/2, CUC2), bHLH 
(basic Helix-Loop-Helix), and WRKY (named after the 
WRKY amino acid motif) [46]. 

In our study, AP2-EREBP, MYB, and GRAS were the 
three major TF families, containing six, five, and five 
genes of the 32 cold-responsive TFs, respectively. The 
AP2-EREBP family plays a major role in the early stages 
of the cold response, as evidenced by many well-charac- 
terized DREBs/CBFs cold regulatory pathway. CBF pro- 
teins, belonging to A-1 subfamily of ERF/AP2 TF family, 
are major regulators that function in activating cold- 
regulated effectors in Arabidopsis and other plants [15]. 
In our study, all members of AP2-EREBP family are 
grouped to different subfamilies of ERF/AP2 TF family. 
For instance, RAP2.11, ERF2, AP2 domain transcription 
factor (CUST_20765), and RAP2A encode members of 
B-6, B-3, A-5, and A-6 subfamilies of ERF/AP2 TF 
family, respectively. Some members have been reported 
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Table 2 Responsive transcription factors during the process of low temperature treatment in cassava apical shoots 



Family 


Probe Name 


Description 


4 h/0 h 


9 h/0 h 
log2 Ratio 


9 h/4 h 


AGI Locus 


E Value 


AP2-EREBP 


CUST_ 18605 


RAP2.1 1 


2.96 


4.59 




AT5G 19790 


1 .40E-06 




CUST_2332 


AP2 domain transcription factor 




4.02 




AT4G34410 


1 .60E-23 




CUST 14955 


ERF2 (ethylene response transcription activator) 






3.35 


AT5G47220 


1.50E-33 




CUST 20765 


AP2 domain transcription factor 




2.02 




AT1G21910 


4.70E-29 




CUST_5161 


RAP2.4 






1.15 


ATI G78080 


8.30E-32 




CUST 2192 


ERF9 transcription repressor 


-1.51 






AT5G44210 


l.OOE-21 


AUX/IAA 


CUST 13515 


IAA16 transcription factor 




-1.55 


-1.10 


AT3G04730 


4.30E-77 


ABI3-VP1 


CUST_3699 


B3 family transcriptional factor 






1.27 


AT4G01580 


1.80E-17 


bZIP 


CUST_1 0706 


ABA-responsive element binding protein 3 




-1.18 




AT3G56850 


1.20E-35 




CUST_1 8065 


ABA-responsive element binding protein 3 




-1.22 




AT3G56850 


2.30E-20 


C2H2 


CUST_4530 


ZFP5 (ZINC FINGER PROTEIN 5) 




1.50 




ATI G 10480 


9.70E-19 




CUST 12461 


C2H2-type zinc finger ZAT-12 like 






2.41 


AT2G28710 


1 .60E-25 


C3H 


CUST_18218 


zinc finger (CCCH-type) family protein 




-1.05 




AT4G29190 


5.00E-77 


G2-like 


CUST_13281 


PCLl (PH\TOCLOCK 1) 




1.62 




AT3G46640 


1.60E-21 




CUST_3593 


PCLl (PH\TOCLOCK 1) 




1.55 




AT3G46640 


4.90E-35 


GRAS 


CUST 13235 


GAI transcription factor 




-1.07 




ATI G 14920 


2.00E-10^ 




CUST 7817 


scarecrow transcription factor 




-1.09 




AT5G66770 


1 .50E-39 




CUST_12751 


SCL5 transcription factor 




-1.39 




AT1G50600 


2.90E-10 




CUST_10812 


scarecrow transcription factor 




-1.52 




AT5G66770 


8.90E-38 




CUST_6129 


scarecrow-like transcription factor 9 






-1.00 


AT2G37650 


2.90E-15 


HSF 


CUST 18574 


heat shock transcription factor A8 




2.14 


1.74 


AT1G67970 


1.20E-17 




CUST 14028 


heat shock transcription factor CI 




1.35 


1.54 


AT3G24520 


9.90E-77 


JUMONJI 


CUST 15293 


jumonji domain transcription factor 




1.85 




AT4G00990 


3.50E-09 


LIM 


CUST 14929 


WLIMl transcription factor 




1.21 




ATI G 10200 


2.20E-86 


MYB 


CUST_9065 


MYB73 






2.01 


AT4G37260 


9.60E-25 




CUST_6759 


myb family transcription factor 




-1.49 




AT1G49010 


9.30E-41 


MYB-related 


CUST_10013 


CPC (CAPRICE) transcription factor 




-1.68 


-1.04 


AT2G46410 


4.10E-12 


NAC 


CUST_1 7347 


NAC domain transcription factor 




-1.34 




AT5G13180 


2.60E-14 


Sigma70-like 


CUST_14809 


SIGMA FACTOR B 




-1.03 




AT1G08540 


5.30E-20 


TCP 


CUST_7714 


PTFl (PLASTID TRANSCRIPTION FACTOR 1) 




-1.26 




AT3G02150 


2.00E-12 


WRKY 


CUST_16901 


WRKY7 transcription factor 




-1.50 




AT4G24240 


2.60E-34 


Unclassified 


CUST_19860 


salt tolerance during germination 1 


-1.38 


-1.17 




AT4G31720 


5.10E-09 



Notes: Positive and negative values of Log2 Ratio are either up- or down-regulated genes In the three pair-comparisons. No significant fold changes are 
Indicated as 



to show response to cold stress in Arabidopsis, such as 
RAP2.11 and RAP2A [47,48]. However, others were 
firstly identified to be involved in cold response, such as 
AP2 domain TF (CUST_2332) and ERF2, Importantly, 
most members of this family were highly induced at 9 h 
during cold treatment, which was confirmed by real- 
time RT-PCR (Figure 4A), suggesting a relatively delay 
in triggering transcriptional cascades in the cold 
response of cassava. Only the ERF9 transcription repres- 
sor was down-regulated after 4 h under cold stress. An 
earlier study had reported that its homologous gene 
(AT5G44210) of ERF9 also experienced a significant 
decrease in its expression after exposure to heat and 
freezing in Arabidopsis [49]. One of CBF homologous 
gene (CUST_6694, hereafter C5F-like) was found to be 



up-regulated following 9 h cold treatment. Due to a 
high P value (0.99), it was excluded from the analysis 
under our standard {P value < 0.05). Such scenario indi- 
cates that the timely response of DREB/CBF pathway to 
cold might be variable among different species. Indeed, 
its cold-induced expression at 9 h following cold treat- 
ment, together with most of the AP2-EREBP family fac- 
tors, was verified by real-time RT-PCR (Figure 4A, B). 
Therefore, based on our analysis, we speculated that 
cassava might mount a relatively slow mobilization of 
stress-related regulators and their downstream genes 
after exposure to cold, rendering cassava vulnerable to 
cold stress. Additionally, the transcript of an Inducer of 
CBF Expression 1 (ICE 1) -like gene of cassava was found 
to be unchanged (Figure 4B). In Arabidopsis and other 
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Figure 4 Expression patterns of important transcripts in response to cold stress at 0, 4, and 9 h by real-time RT-PCR analysis. (A) 

Members of AP2-EREBP transcription factor family. (B) C8F-lil<e gene and its upstream regulons. (C) ProDH, P5CS, and ASPS engaged in 'Proline 
metabolism'. (D) AMYl and TPS participated in 'Starch and sucrose metabolism'. The double asterisks indicate a statistically significant difference 
(p < 0.01) for the data of the stress-treated samples compared to those of the unstressed samples. 



plants, ICEs have been reported to be constitutively 
expressed [50]. 

A comprehensive expression analysis of MYB TFs 
demonstrated that almost all MYB TFs are responsive 
to stresses or hormones [51]. In our study, five MYB or 
MYB-related family members showed differential 
expression, with 3 cold-inducible and 2 cold-repressed. 
Two transcripts, encoding PCLl with a single myb 
DNA-binding domain, were required for circadian 
rhythms [52], and the other three MYB/MYB-related 
transcription factors were responsive to a diverse of hor- 
mone treatments in Arabidopsis [53]. Differential 
expression of MYB TFs implies that other environmen- 
tal or hormonal pathways may be involved in cold 
response in cassava. Furthermore, the up-regulation of 
the CUST_3593 transcript {PCLl) was also validated by 
real-time RT-PCR. Members of the GRAS gene family 
encode transcriptional regulators that have diverse func- 
tions in plant growth and development, such as gibber- 
ellin signal transduction, root radial patterning, and 



axillary meristem formation [54]. Although the Arabi- 
dopsis genome encodes at least 33 GRAS protein family 
members, few GRAS proteins have been characterized 
thus far [55]. It is worth noting that all GRAS family 
members were down-regulated on the array (Table 2), 
suggesting that cold stress inhibit plant growth and con- 
serve energy to adapt to the adverse environment [56]. 
Similarly, HSF TFs were reported to be involved in the 
response to heat and other abiotic stresses [38,57]. The 
microarray analysis also showed the inducible expression 
of heat shock transcription factor A8 and CI following 
cold treatment (Table 2). Interestingly, AS was also 
found to be induced in response to low temperatures in 
potato and Arabidopsis [36]. The induction of heat 
shock TFs and heat shock proteins (Table 1) revealed 
that cold and heat stresses may share common respon- 
sive elements. The majority of WRKY family TFs is 
known to be responsive to biotic and/or abiotic stress, 
despite the fact that most research has focused on the 
role of these genes in plant-pathogen interactions [58]. 
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WRKY7, known to encode a novel CaM-binding tran- 
scription factor of the WRKY group Ild in Arabidopsis 
[59], was also found to have differential expression in 
our study. This might indicate a crosstalk between cold 
stress response and plant-pathogen interaction in 
cassava. 

In addition to the above-mentioned TFs, we found 
two AREBs that were cold-repressed in our study. These 
factors belonged to group A bZIP transcription factors, 
which play a role in plant pathogen responses, light sig- 
naling, and ABA and abiotic stress signaling, and were 
induced by salt and ABA in Arabidopsis [60]. The zinc 
finger family has been demonstrated to be involved in 
the cold response in Arabidopsis [38]. In our study, we 
found 3 probes that fell into 2 subfamilies of the zinc 
finger family (C2H2 and C3H), including 2 up-regulated 
and 1 down-regulated genes. Zinc finger proteins 
involve in ROS and abiotic stress signaling in Arabidop- 
sis [61]. The C2H2-type zinc finger protein gene 
ZAT12'\ike was induced not only by cold stress, but 
also by PEG and salt stress, as illustrated by our real- 
time RT-PCR results (Figure 4B). This indicates that 
cassava needs to trigger the expression of ROS-scaven- 
ging genes to adapt oxidative stress. The plant-specific 
NAC transcription factor family has been implicated in 
plant development processes, such as shoot apical mer- 
istem (SAM) maintenance and organ differentiation 
[62], as well as biotic and abiotic stress responses [63]. 
We found only one transcript encoding for the NAC 
domain transcription factor, which was down-regulated 
on the array. In addition, five novel transcription factor 
families (ABI3-VP1, JUMONJI, LIM, Sigma70-like, and 
TCP) were also identified. Their homologous genes in 
other plant species have not yet been reported in 
response to cold stress, suggesting that these genes 
might be specific to cassava and are attractive targets 
for further functional characterization. 
Cold-responsive genes related to 'Chloroplasf 
The category of 'Chloroplast' was an abundant and 
important GO Slim. It accounted for 21.06% of total 
GOs in 'Cellular Component' (Figure 3C, Additional file 
11), indicating that a great number of chloroplast-asso- 
ciated genes changed their expression to adapt to cold. 
This result suggests that the chloroplast might be one of 
the major organelles affected, as supported by the obser- 
vation of damaged chloroplasts in cold-stressed leaves 
(Figure IC). Although the differential expression of 
chloroplast-related genes might also be affected by other 
environmental changes on cassava plants during the 
treatment, i.e. light change, low temperature was the 
dominant factor in this study. 

Chloroplast is an important organelle unique to plant 
cells and is the site of photosynthesis. In previous reports, 
cold stress could reduce the rate of photosynthesis by 



interfering with the function of many photosynthesis- 
related proteins [64]. In our study, 43 out of 319 cold- 
responsive transcripts were assigned to the 'Chloroplast' 
class, including 22 up-regulated and 21 down-regulated 
genes (Table 3). Among them, several genes encoding 
protein localized to thylakoid were identified. For exam- 
ples, CABl (chlorophyll a/b binding protein 1), a lipid- 
associated family protein, a thylakoid leumenal 20 kDa 
protein, and a hypothetical protein (CUST_7076) were all 
up-regulated under cold stress. In contrast, LHCA4 
(photosystem I light harvesting complex gene 4), PsaL 
(photosystem I subunit L), and the rubredoxin family 
protein were dramatically down-regulated. Two of them 
{LHCA4 and PsaL) have also been proposed to be 
involved in photosynthesis in Arabidopsis. These results 
suggested that differential expression of chloroplast 
genes, especially those encoding for thylakoid associated 
proteins, by cold stress might lead to the processes of 
thylakoid distortion, chloroplast malfunction, and photo- 
synthesis inhibition. 

KEGG pathway analysis of cold-responsive genes 

To determine whether the cold stress responsive genes 
engaged in specific pathways, 291 Arabidopsis AGI loci 
representing the 319 DEGs were used as objects to 
search against KEGG pathway maps in Arabidopsis 
thaliana. Finally, 44 related pathways were identified 
(Figure 5, Additional file 12). Several interesting and 
important pathways, including 'Plant hormone signal 
transduction' (ath04075), 'Plant-pathogen interaction' 
(ath04626), 'Phenylalanine metabolism' (ath00360), 
'Fatty acid biosynthesis' (ath00061), and 'Starch and 
sucrose metabolism' (athOOSOO), were involved and func- 
tion in the early response to cold stress. 

'Plant hormone signal transduction' (ath04075) com- 
prised 8 genes of the 44 cold-regulated and pathway-hit 
genes on our array. In this pathway, the transcripts of 
GAl, SAUR (auxin-responsive protein), EBFl (EIN3- 
binding F box protein 1), EIN4 (ethylene insensitive 4), 
and AREBS were determined to be down-regulated, and 
the transcript of ERF2 was identified as being up-regu- 
lated. Their homologous genes have been reported to be 
implicated in GA, auxin, ethylene (ETH), and ABA- 
mediate hormone signal transduction, respectively 
[65,66]. In plants, these hormones play crucial roles in a 
diverse set of developmental processes, as well as biotic 
and abiotic stresses [67]. For example, ABA accumulates 
in response to abiotic stresses, such as cold, salt, and 
drought [68]. Among these hormone-related genes in 
Arabidopsis, GAI has been also reported to decrease its 
transcript levels in response to cold [38], and ERF2 has 
been identified in the response to chitin, a plant-defense 
elicitor [69]. Previous studies have suggested that the 
differential regulation of genes involved in hormone 
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Table 3 Low temperature responsive transcripts related to 'Chloroplast' 



Probe Name 


Description 


4 h/0 h 


9 h/0 h 
Log2 Ratio 


9 h/4 h 


AGI Locus 


E Value 


CUST 12882 


CCL (CCR-LIKE) 




4.19 


3.26 


AT3G26740 


1 .90E-25 


CUST_6863 


FABl (fatty acid biosynthesis 1) 


3.45 


4.17 




ATI G74960 


2.10E-13 


CUST_9556 


CABl (chloropliyll a/b binding protein 1) 


2.84 


4.06 


1.22 


ATI G29930 


9.00E-134 


CUST 141 18 


SSI2 (fatty acid biosynthesis 2) 


2.31 


3.69 


1.38 


AT2G43710 


1.80E-57 


CUST_16917 


lipid-associated family protein 


2.09 


3.48 




AT4G39730 


5.50E-44 


CUST_222 


PDSl (PHnOENE DESATURATION 1) 




2.47 


1.74 


AT1G06570 


2.00E-13 


CUST 19871 


dehydrodolichyl diphosphate synthase 




2.21 




AT5G58770 


1 .60E-42 


CUST 3426 


DNA binding 




2.15 


1.40 


AT2G37020 


1 .60E-44 


CUST_1 1 923 


MEE14 (maternal effect embryo arrest 14) 




1.85 


1.81 


AT2G 15890 


3.40E-42 


CUST 8767 


hydrolase, alpha/beta fold family protein 




1.76 




AT5G 13800 


2.70E-33 


CUST_4530 


ZFP5 (ZINC FINGER PROTEIN 5) 




1.50 




ATI G 10480 


9.70E-19 


CUST 15215 


thylakoid lumenal 20 kDa protein 




1.45 


1.13 


AT3G56650 


2.50E-22 


CUST_1 0800 


hydrolase 


1.03 


1.38 




AT2G35450 


4.30E-07 


CUST_6415 


pentatricopeptide (PPR) repeat-containing protein 




1.28 




AT3G02650 


9.00E-20 


CUST 9148 


AMY3 (alpha-amylase-like 3) 




1.24 




ATI G69830 


2.80E-53 


CUST 20805 


BIGYIN 


1.36 


1.22 




AT3G57090 


3.40E-23 


CUST 3083 


hypothetical protein 




1.13 




AT5G59400 


l.lOE-46 


CUST_3051 


CAT2 (catalase2) 




1.1 1 


1.01 


AT4G35090 


7.60E-76 


CUST_4761 


hypothetical protein 




1.07 




AT4G 10000 


6.40E-10 


CUST_20137 


NFU3 (NFU domain protein 3) 




1.04 




AT4G25910 


4.30E-17 


CUST_7076 


hypothetical protein 




1.01 




AT2G42130 


1 .30E-40 


CUST_4922 


CHUPl (CHLOROPLAST UNUSUAL POSITIONING 1) 




-1.00 




AT3G25690 


1 .60E-09 


CUST_ 14809 


SIGB (SIGMA FACTOR B) 




-1.03 




ATI G08540 


5.30E-20 


CUST_14166 


GSTU19 (glutathione transferase) 




-1.12 




ATI G78380 


1 .90E-70 


CUST_1453 


CHUPl (CHLOROPLAST UNUSUAL POSITIONING 1) 




-1.13 




AT3G25690 


2.60E-95 


CUST_12230 


amine oxidase family protein 




-1.22 




AT3G09580 


4.10E-90 


CUST_7714 


PTFl (PLASTID TRANSCRIPTION FACTOR 1) 




-1.26 




AT3G02150 


2.00E-12 


CUST 18664 


rubredoxin family protein 




-1.28 




ATI G54500 


2.90E-43 


CUST 1805 


hypothetical protein 




-1.29 




AT5G40470 


3.90E-15 


CUST_16156 


dehydrodolichyl diphosphate synthase 




-1.29 




AT5G58770 


1.20E-14 


CUST 15952 


CHUPl (chloroplast unusual positioning 1) 




-1.37 




AT3G25690 


1 .70E-69 


CUST_1 9234 


LHCA4 (photosystem 1 light harvesting complex gene 4) 




-1.39 




AT3G47470 


8.80E-51 


CUST_2381 


NCED4 (nine-cis-epoxycarotenoid dioxygenase 4) 




-1.45 




AT4G19170 


1 .OOE-63 


CUST_15061 


hypothetical protein 




-1.48 


-1.08 


AT3G22970 


1 .80E-44 


CUST_5552 


pentatricopeptide (PPR) repeat-containing protein 




-1.49 




AT3G42630 


1 .40E-63 


CUST_10020 


thioredoxin family protein 




-1.52 




AT1G08570 


1 .40E-43 


CUST_2294 


ATOFP15/OFP15 




-1.55 


-1.09 


AT2G36050 


1.40E-28 


CUST_7622 


DNAJ heat shock protein, putative 




-1.57 




AT2G 17880 


9.20E-25 


CUST_5765 


DVL13/RTFL2 (ROTUNDIFOLIA LIKE 2) 




-1.62 




AT2G29125 


2.50E-16 


CUST_3996 


hypothetical protein 




-1.64 


-1.05 


AT1G68430 


4.50E-13 


CUST_8087 


xylulose kinase, putative 




-1.98 




AT2G21370 


1 .60E-56 


CUST_6408 


PSAL (photosystem 1 subunit L) 


-2.47 


-2.31 




AT4G 12800 


9.80E-48 



Notes: Positive and negative values of Log2 Ratio are either up- or down-regulated genes In the three pair-comparisons. No significant fold changes are 
Indicated as 



signal transduction might play key roles in the early cold 
stress response [38]. 

Two other important pathways, including 'Phenylala- 
nine metabolism' (ath00360) and 'Plant-pathogen inter- 
action' (ath04626), were also found in our study to be 
regulated. PDSl (phytoene desaturation 1), ASPS 



(Aspartate aminotransferase 3), and peroxidase 72, asso- 
ciated with 'phenylalanine metabolism', were all highly 
accumulated in response to cold. As noted in previous 
studies, PDSl has been reported to be engaged in the 
carotenoid and plastoquinone biosynthetic process [70]. 
ASPS, also participating in proline metabolism, has been 
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others (28) 
44% 
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3% 
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Arginine and 

proline 
metabolism 
(2) 
3% 
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Photosynthesis (2) 

3% 
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3% 




ath04075 Plant hormone signal 
transduction (6) 

ath04626 Plant- 
pathogen 
-interaction (3) 

5% 



[h00360 Phenylalanine 
metabolism (3) 
5% 



ath00061 Fatty acid 
biosynthesis (2) 
3% 



ath00270 Cysteine and 
methionine metabolism (2) 

^''o ath03022 Basal 
transcription factors (2) 
3% 



ath00380 Tryptophan 
metabolism (2) 

3% 



ath00250 Alanine, aspartate and 
glutamate metabolism (2) 
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sucrose metabolism (2) 

3% 



ath00350 Tyrosine metabolism (2) 

3% 



Figure 5 KEGG pathway maps of cold-responsive genes. A total of 44 pathways were identified tlirougli KEGG mapping. Different colors 
represent the pathway entries and pathway names. The number of genes represented by AGI loci involved in each pathway is labeled in the 
parentheses. 



shown to be associated with leaf senescence and to be 
induced in response to darkness and ethylene-induced 
artificial senescence [71]. In Arabidopsis and other plants, 
proline levels are mainly determined by balance of bio- 
synthetic and catabolic pathways, controlled by P5CS 
(^^-pyrroline-5-carboxylate synthetase) and ProDH (pro- 
line dehydrogenase) genes, respectively [72]. Up-regula- 
tion of cassava homologous genes PSCS and ASPS, and 
down-regulation of ProDH, which were confirmed by 
real-time RT-PCR (Figure 4C), might correlate with the 
high accumulation of proline (Figure IE). Peroxidase 72, 
a member of the Class III peroxidases, is involved in the 
response to oxidation stress and potassium resupply [73]. 
Additionally, transcripts from MAP Kinase 4, UNE14 
(unfertilized embryo sac 14), and TCH2, members of the 
'Plant-pathogen interaction' pathway, were also induced 
by cold (Table 4). MAP Kinase 4 and TCH2 have been 
shown to be involved in response to several stresses, 
including cold, salinity, and heat in Arabidopsis [39]; 
whereas UNE14 also encodes a calcium ion binding pro- 
tein. These results indicate that calcium ion binding pro- 
teins and kinase-mediated signal transduction may play 
pivotal roles in cold stress response in cassava. 

In addition, 'Fatty acid biosynthesis' (ath00061) and 
'Starch and sucrose metabolism' (athOOSOO) were also 
identified through the pathway analysis. FABl (fatty acid 
biosynthesis 1) and SSI2, which were also assigned to 
'Chloroplast' category, were involved in 'Fatty acid bio- 
synthesis'. Earlier studies in Arabidopsis indicated that 
FABl is correlated with chilling stress [74], whereas 
SSI2 has been shown to respond to biotic attacks. 



including viruses, insects, and bacteria [75]. Both were 
strongly induced during the entire cold treatment pro- 
cess, indicating that membrane modification may occur 
during the early cold stress response. In Arabidopsis, 
TPS7 [alpha, alpha-trehalose-phosphate synthase (UDP- 
forming) /transferase] has been reported to participate in 
the trehalose biosynthetic process. AMYl (alpha-amy- 
lase-like 1) has also been shown to take part in the 
degradation of starch to sugar and to be induced by bio- 
tic and abiotic stress [76]. In our study, TPS7 and AMYl 
were identified as being involved in 'Starch and sucrose 
metabolism'. Their cold-induced expression was con- 
firmed by real-time RT-PCR analysis (Figure 4D). In 
cold-stressed leaves, starch granules were found to have 
disappeared (Figure IC). Even though their total sugar 
content remained unchanged after 4 h and was lower 
after 9 h of treatment, an obvious increase was detected 
after 24 h (Figure 6A). Interestingly, sucrose was signifi- 
cantly accumulated and reached the highest level in 9 h; 
whereas glucose levels were only slightly increased 
under cold stress (Figure 6B). These results suggested 
that soluble sugars might act as both osmolytes and sig- 
nal molecules in the cold response of cassava, similar to 
other plant species [77]. It also provided insight into the 
hypothesis that starch degradation provides conserve 
energy for plants under adverse conditions. Such transi- 
tion from synthesis metabolism to degradation metabo- 
lism would also explain the reduced productivity 
observed in plants under stress conditions. 

In short, a variety of interesting and important path- 
ways are involved and function in the cellular response 
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Table 4 Validation of selected microarray-based gene expression by real-time RT-PCR analysis 



Probe Name 


Description 


Microarray 




qRT-PCR 




AGI Locus 


Category 






4 h/0 h 


9 h/0 h 


2 h/0 h 


4 h/0 h 


9 h/0 h 


24 h/0 h 






CUST_9007 


JAZ7 (jasmonate-ZIM-domain protein 7) 






2.34 


1 Qf\ 


^ 7^ 
J./ J 


2.82 


AT2G34600 


UK 


CUST_5663 


Unknown protein 


1 f^Q 
1 .Oy 




1.20 


1 

1 .DO 


9 A^ 
Z.Oj 


4.19 






CUST_5990 


phosplioric monoester hydrolase 




Q 77 
J.I 1 


0.54 


n 1 1 
U. 1 1 




2.19 


ATI Gl 7710 


UK 


CUST 4383 


gibberellin 2-beta-dioxygenase 




J.DO 


5.80 


O.OO 


A CM 


4.62 


ATI G30040 


RS 


CUST 13673 


protein transmembrane transporter 




j.D 1 


1.41 


u.zu 


1 71 
1 ./ 1 


2.08 


AT4G16160 


UK 


CUST 16591 


Unl<nown protein 


z.zo 


J.J J 


1.41 


1 .DZ 


1 ^ZL 
1 .D'-t 


2.61 


AT5G37770 


RS 


CUST_1 1786 


ubiquitin-protein ligase 






1.14 


1 1 Q 
1.13 


1 

1 .J J 


1.91 


AT1G68050 


RS, CR 


CUST_4231 


sodium/calcium exclianger protein 




Q 1 1 
J. \ \ 


0.35 


n AH 

U.OU 


1 Q7 
1 .0/ 


1.32 


AT1G53210 


UK 


CUST 20765 


AP2 domain transcription factor 




2.02 


0.23 


0.42 


2.41 


1.70 


AT1G21910 


TF 


CUST_12461 


C2H2-type zinc finger ZAT-12 lil<e 




1.97 


1.23 


1.39 


1.10 


2.57 


AT2G28710 


TF 


CUST_19994 


protein pliospliatase 2C 




1.71 


0.68 


0.39 


1.12 


0.28 


AT3G 12620 


PM 


CUST_3593 


PCLl (PH\TOCLOCK 1) 




1.55 


1.17 


0.95 


1.96 


2.40 


AT3G46640 


TF 


CUST_17985 


TCH2 (calcium ion binding protein) 




1.53 


2.66 


0.59 


1.41 


1.79 


AT5G37770 


PP 


CUST_14166 


GSTU19 (glutathione transferase) 




-1.12 


-0.82 


-0.67 


-1.05 


-1.00 


AT1G78380 


RS, CH 


CUST_1923 


S-adenosylmethionine decarboxylase 




-1.85 


-0.91 


-1.56 


-1.89 


-2.63 


AT3G02470 


AP 


CUST_7527 


Unknown protein 


-1.47 


-1.85 


-0.43 


-1.36 


0.92 


-1.06 




UK 


CUST_19213 


Unknown protein 


-1.64 


-2.59 


-0.77 


-1.27 


-1.28 


-1.39 




UK 


CUST_5741 


hypothetical protein 


-2.61 


-3.05 


-1.94 


-2.39 


-3.04 


-3.48 


ATI G 13360 


UK 



Notes: All data are shown log2 Ratio, and positive and negative values of Log2 Ratio are either up- or down-regulated genes In the three pair-comparisons. No 
significant fold changes are Indicated as Values of real-time RT-PCR analysis showing In bold font Indicate Inconsistent data compared with microarray 
analysis. AP: 'Arglnlne and proline metablism', CH: 'Chloroplast', CR: 'CIrcadlan rhythm-plant', PM: 'Protein metabolism', PP: 'Plant-pathogen Interaction', RS: 
'Response to abiotic and biotic stimulus' and 'Response to stress', TF: 'Transcription factor activity', UK: Unknown. 



to cold stress. Despite many pathways, e.g., 'Plant hor- 
mone signal transduction' (ath04075), 'Fatty acid bio- 
synthesis' (ath00061), and 'Starch and sucrose 
metabolism' (athOOSOO), have been well-demonstrated in 
other plant species, regulation of different pathways or 
gene has been observed, illustrating that a complex and 
specific network is involved in the early cold response in 
cassava. 

Validation of microarray results by real-time RT-PCR 

To validate our microarray data, real-time RT-PCR ana- 
lysis was performed on 18 selected differentially 
expressed transcripts, including 13 up-regulated and 5 
down-regulated genes (Table 4). These genes belong to 
divergent functional categories or pathways. For 
instance, four genes [gibberellin 2-beta-dioxygenase, 
unknown protein (CUST_16591), ubiquitin-protein 
ligase, and GSTU19\ were involved in 'Response to abio- 
tic and biotic stimulus' and 'Response to stress'. One 
gene (protein phosphatase 2C) was implicated in 'Pro- 
tein metabolism,' and other three genes [AP2 domain 
transcription factor (CUST_20765), C2H2-type zinc fin- 
ger (ZAr72-like), and PCLl (CUST_3593)] were 
grouped as transcription factors. In addition, 3 genes 
(ubiquitin-protein ligase, TCH2, and S-adenosylmethio- 
nine decarboxylase) were assigned to 'Circadian rhythm- 
plant' (ath04712), 'Plant-pathogen interaction', and 
'Arginine and proline metabolism' (ath00330). 



respectively. Seven genes that are unclassified in the GO 
Slims or pathways but displayed significant changes in 
transcript abundance were also tested. For these genes, 
the fold change (-^^Ct) measured by real-time RT-PCR 
and by microarray was highly consistent (72.2% and 
94.4% for 4 h/0 h and 9 h/0 h, respectively). Their 
expression kinetics from the real-time RT-PCR results 
was similar to those of the microarray analysis. These 
results re-confirmed the accuracy of our microarray 
data. 

Genes responsive to different stresses and their tissue- 
specific expressions 

The cold stress-signaling pathway may interact with 
other signaling systems of, for example, ABA, salt, and 
drought [19]. The expression patterns of the 13 cold- 
responsive genes were analyzed by real-time RT-PCR 
using cassava in vitro shoot cultures treated with 100 
(iM ABA, 25% PEG, and 250 mM NaCl for 6 h. Genes 
with fold changes larger than 4 were defined as strongly 
responsive to the stresses. Several cold-responsive genes 
were also strongly induced by ABA, PEG, or salt stress 
(e.g., CUST_16591, novel protein with unknown func- 
tion; CUST_12461, C2H2-type zinc finger, ZATi2-like; 
CUST_4383, gibberellin 2-beta-dioxygenase) (Table 5). 
These results indicated that dynamic crosstalk exists 
among signaling pathways related to cold, drought, and 
salt in cassava because all of these stresses cause cellular 
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Sucrose Glucose 

Figure 6 Contents of total soluble sugars (A) and sucrose and 
glucose (B) in cassava leaves subjected to cold (7°C) for 0, 4, 9, 

and 24 h. The single and double asterisks indicate a statistically 
significant difference at p < 0.05 and p < 0.01, respectively, for the 
data of the stress-treated samples compared to those of the 
unstressed samples. The mean values are calculated from three 
biological replicates; the error bars represent the standard error of 
the mean (SEM). 



dehydration [20]. Genes uniquely responsive to cold 
stress drew more attention in our study, such as TCH2 
(CUST_17985), AP2 domain transcription factor 
(CUST_20765), and JAZ7 (jasmonate-ZIM-domain pro- 
tein 7, CUST_9007) (Table 5). Both TCH2 and JAZ7 
were involved in the 'Plant-pathogen interaction' in Ara- 
bidopsis (Additional file 12) and rice, suggesting that 
chilling-induced processes share some common features 
with the defense mechanism against pathogens in plants 
[58]. In addition, protein transmembrane transporter 
(CUST_13673), which was strongly induced by ABA and 
salt stress, was also believed to function in early signal 
reception. 

To further analyze the expression patterns of these 
stress-regulated genes in different tissues, real-time RT- 
PCR analysis was conducted. mRNAs were extracted 
from the apical buds (AP), fibrous roots (FR), mature 
leaves (ML), stem cambia (SC), young leaves (YL), and 
young stems (YS) of 3-month-old cassava plants grown 
in a greenhouse. Interestingly, 7 tested genes were pri- 
marily expressed in FR, ML, SC, and YS rather than in 
AP and YL (Table 6). This tissue-specific expression of 
cold-responsive genes might explain why the apical 
shoots, including AP and YL, were more susceptible to 
cold stress than ML, but this difference remains to be 
further determined. 

Changes in H2O2 content and ROS scavenging enzyme 
activities 

As previous reported, chilling conditions may lead to an 
accumulation of ROS such as hydrogen peroxide (H2O2) 
and superoxide radical (O2') [78]. To analyze and visua- 
lize H2O2 produced in the cassava leaves subjected to 
cold stress (7°C), diaminobenzidine (DAB) was used to 



Table 5 Real-time RT-PCR verification of cold responsive gene under the different abiotic stresses and ABA treatment 

Probe Name Description ABA PEG NaCI Cold AG! Locus E Value 



Log2 Ratio 



CUST_12461 


C2H2-type zinc finger ZAT-12 like 


1.45 


2.26 


2.63 


3.89 


AT2G28710 


1 .60E-25 


CUST_16591 


Unl<nown protein 


2.52 


4.91 


2.49 


2.30 






CUST_4383 


gibberellin 2-beta-dioxygenase 


1.12 


1.71 


2.52 


2.65 


AT1G30040 


4.70E-51 


CUST_17985 


TCH2 (calcium ion binding protein) 


-0.72 


-0.62 


-0.41 


2.75 


AT5G37770 


2.10E-39 


CUST_20765 


AP2 domain transcription factor 


1.08 


-0.61 


0.38 


2.16 


AT1G21910 


4.70E-29 


CUST_9007 


JAZ7 (jasmonate-ZIM-domain protein 7) 


0.11 


-0.57 


-0.43 


2.16 


AT2G34600 


2.70E-10 


CUST_9556 


CABl (chlorophyll a/b binding protein 1) 


1.24 


0.21 


0.65 


1.94 


AT1G29930 


9.00E-134 


CUST_13673 


protein transmembrane transporter 


3.22 


1.85 


4.04 


1.73 


AT4G16160 


3.40E-55 


CUST_5990 


phosphoric monoester hydrolase 


1.21 


-1.58 


0.09 


1.64 


ATI G1 7710 


8.60E-90 


CUST_14166 


GSTU19 (glutathione transferase) 


1.00 


1.66 


1.01 


-0.74 


AT1G78380 


1 .90E-70 


CUST_19213 


Unknown protein 


0.41 


0.00 


-0.13 


-1.43 






CUST_7527 


Unknown protein 


-0.11 


-1.21 


-1.83 


-0.70 






CUST_1923 


S-adenosylmethionine decarboxylase 


0.59 


-0.90 


-0.64 


-0.70 


AT3G02470 


1 .40E-08 



Notes: Thirteen selected cold responsive genes were analyzed under different treatments, Including 100 jjM ABA, 25% PEG, 200 mM NaCI, and cold stress (7°C) 
for 6 h. Values that were dramatically higher (Log2 Ratio > 2) than those of control are shown In bold font. 
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Table 6 Expression patterns of stress responsive genes in different tissues and organs 



Probe Name 


Description 


AP 


FR 


ML 


sc 


YL 


YS 








Log2 Ratio 






CUST_9007 


JAZ7 (jasmonate-ZIM-domain protein 7) 


1.00 


4.03 


4.29 


1.48 


1.01 


8.05 


CUST_12461 


C2H2-type zinc finger ZAT-12 lil<e 


1.00 


6.76 


6.72 


4.02 


1.93 


5.17 


CUST_5990 


pliosplioric monoester liydrolase 


1.00 


3.32 


8.24 


5.32 


0.23 


3.82 


CUST_17985 


TCH2 (calcium ion binding protein) 


1.00 


3.24 


3.24 


2.77 


1.44 


4.57 


CUST_4383 


gibberellin 2-beta-dioxygenase 


1.00 


1.28 


1.40 


0.88 


0.38 


4.78 


CUST_13673 


protein transmembrane transporter 


1.00 


6.16 


3.96 


2.93 


-0.04 


1.77 


CUST_16591 


Unl<nown protein 


1.00 


2.55 


0.74 


3.09 


-0.21 


-0.34 



Notes: AP, apical buds; FR, fibrous roots; ML, mature leaves; SC, stem cambia; YL, young leaves; YS, young stems. All expression data shown were compared with 
AP. Values that were dramatically higher (Log2 Ratio > 2) than those of control are shown In bold font. 



Stain the first or second fully expanded leaves of in vitro 
specimens and greenhouse-grown cassava plants. In 
unstressed leaves, few DAB polymerizations could be 
observed; in contrast, there was a significant increase in 
DAB polymerization detected as early as 4 h following 
cold treatment that continued up to 24 h (Figure 7A). 
This result was consistent with the quantification of 
H2O2 content (Figure 7B), indicating that the oxidative 
stress may exert a toxic effect on cassava to adapt or 
survive under cold stress conditions. 

Plants, as well as other organisms, have evolved anti- 
oxidant systems to protect themselves against toxic spe- 
cies of oxygen. ROS scavenging enzymes, including 
catalase (CAT), superoxide dismutase (SOD), and glu- 
tathione transferase (GST), have been demonstrated to 
play key roles in the removal of ROS. In our DEGs, sev- 
eral genes encoding ROS scavenging enzymes were 
recorded to be target to chloroplast, including CAT2 
and GSTU19 (Table 3). In addition, peroxidase 72, and 
other three GSTs {GSTU7, GSTill, and microsomal 
GST) were also identified (Additional file 4). The up- 
regulation of four ROS-scavenging transcripts (CAT2, 
peroxidase 72, GSTl and microsomal GST) was verified 
by real-time RT-PCR (Figure 7C). In Arabidopsis, CAT2 
was reported to be induced by cold stress and GSTU19 
was responsive to many biotic and abiotic stresses, such 
as water deprivation and oxidative stress [79]. To further 
validate the function of ROS scavenging enzymes after 
exposure to cold stress, the activities of CAT and SOD 
in cassava leaves subjected to 7°C were measured. Con- 
sistent with the microarray data, an increase in CAT 
activity was observed after 9 h of cold treatment (Figure 
7D). Similarly, a significant elevation in SOD activity 
was also observed after 4 h of cold stress, and such 
increases continued for 24 h after the initiation of stress 
treatment (Figure 7E). These results including induction 
of ZAT12'\ike gene (Figure 4B) suggested that ROS sig- 
naling pathways, including ROS scavenging enzymes, 
were involved in the ROS detoxification induced by the 
cold stress response in cassava. Similarly, rice and 



maize, as chilling-sensitive tropical crops, could only 
withstand transient and milder cold stress. In japonica 
rice, an oxidative-mediated network has been proposed 
to play a key role in the early response to chilling stress 
and short-term defenses [80], and insufficient antioxi- 
dant defenses are thought to cause maize chilling sensi- 
tivity [81]. Generally, the ultra-structural changes of the 
chloroplast (Figure IC) and oxidative burst were typi- 
cally characteristic of PCD in plants, indicating that 
PCD may be a part of an adaptive mechanism to survive 
stress [82]. 

Conclusions 

Our study presented a genome-wide gene expression 
profiling of cassava subjected to cold stress using the 
microarray technology. Overall, the transcriptomic 
responses of cassava to cold stress were basically consis- 
tent with the changes seen in other plants under abiotic 
stresses [83]. A considerable amount of specific cassava 
genes related to different biological functions were iden- 
tified. For examples, many new stress-responsive genes 
or TFs, such as early flower 4, ERF2, and AP2 domain 
transcription factor (CUST_2332), were found, suggest- 
ing that various regulatory pathways may exist in cas- 
sava together with the well-characterized CBF pathway. 
Based on the comprehensive transcriptomic, real-time 
RT-PCR and physiological analyses, our study supports 
the fact that crosstalks among different signaling sys- 
tems play an important role in regulating the cold stress 
response in tropical plants. Thus, a hypothetical model 
for depicting the components involved in cold response 
networks in cassava could be established (Figure 8). 
Plants may perceive low temperature through cell mem- 
brane receptors, most likely transmembrane proteins, 
and by membrane modification through fatty acid synth- 
esis (e.g., FABl and SSI2). Then, the cold signal trans- 
duction might induce cellular metabolism changes, such 
as the calcium signaling cascade, hormone signal trans- 
duction, and ROS signaling. The calcium signaling path- 
way may stimulate the calcium ion binding protein. 
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0 h 



4h 



9 h 



24 h 





CUST_8448 CUST_3051 CUST_11831 CUST_12421 
(Peroxidase 72) (CAT2) (GSTT1) (microsomal GST) 




Figure 7 Detection of H2O2/ transcript level and enzymatic activities of ROS scavenging genes in cassava subjected to cold (7°C) for 0, 

4, 9, and 24 h. (A) DAB polymerization of in vitro and greenhouse-grown leaves subjected to cold. (B) Quantification of H2O2 content in the 
leaves from the greenhouse-grown plants. (C) The transcript levels of genes encoding for ROS scavenging enzymes in response to cold stress. 
(D, E) Enzyme activities of catalase and superoxide dismutase (SOD) in the leaves of greenhouse-grown plants. The single and double asterisks 
indicate a statistically significant difference at p < 0.05 and p < 0.01, respectively, for the data of the stress-treated samples compared to those 
of the unstressed samples. The mean values are calculated from three biological replicates; the error bars represent the standard error of the 
mean (SEM). Bar = 1 cm. 
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Figure 8 Molecular model of the early cold response in cassava. Two biological processes, namely stress perception and the physiological 
response, are illustrated. After signal reception, stress-activated Ca^^ signaling, ROS signaling, and hormone signaling modulate the expression of 
stress-responsive genes, which include metabolic enzymes, transcription factors, kinases, and ion transporters. The physiological changes that 
manifested as membrane modification, chloroplast malfunction, and starch and sucrose metabolism, as well as amino acid metabolism, cause 
cassava to either have increased cold tolerance or to enter into accelerated programmed cell death. Their balance determines the outcome of 
the stressed cassava plant. Selected up- and down-regulated genes are in red or blue, respectively. The black arrows indicate a positive effect, 
and the black dashed arrows indicate a negative effect. The thick gray arrows show different biological processes in the stress response. AMYl: 
alpha-amylase like 1, AP2-EREBP: APETALA2-Ethylene Responsive Element Binding Protein, AREB3: ABA-responsive element binding protein 3, 
ASP3: Aspartate aminotransferase 3, AUX: auxin, CaM: calmodulin, EBF1: EIN3-binding F box protein 1, EIN4: ethylene insensitive 4, ERF2: ethylene 
response transcription activator, ETH: ethylene, FAB1: fatty acid biosynthesis 1, GA: gibberellins, GAI: gibberellic acid insensitive, GRAS: GAI, RGA, 
and SCR, GSTs: glutathione transferases, HSF: heat shock factor, JA: jasmonate, JAZ7: jasmonate-ZIM-domain protein 7, MARK: Mitogen-activated 
protein kinase, PDSl: Phytoene desaturation 1, SAUR: auxin-responsive protein, SOD: superoxide dismutase, SSI2: fatty acid biosynthesis 2, TPS7: 
trehalose-phosphatase/synthase 7. 



leading to the activation of cascade I<inase activities (e.g., 
MAP kinase 4), which could switch on various cold 
stress-responsive genes and transcription factor family 
proteins (e.g., AP2-EREBP, HSF, and GRAS). Hormone 
(auxin, ABA, GA, ETH, and JA) signal transduction 



could also be activated to alter plant growth status in 
order to adapt to stress condition through the differen- 
tial regulation of their downstream genes (e.g., SAUR, 
AREB3, GAI, ERF2, and JAZ7, respectively). Moreover, 
other metabolisms could also take part in the process. 
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including starch and sucrose metabolism (e.g., TPS7 and 
AMYl), amino acid metabolism (e.g., ASPS and PDSl), 
and the ROS scavenging system (e.g., catalase 2, SOD, 
and GSTs). Apparently, the balance between cell 
damage and tolerance might decide the fate for cold- 
stressed cassava plants (Figure 8). On the other hand, a 
large number of genes (about 37.2% of total 508 DEGs) 
encode proteins of unknown functions. Studying these 
genes may reveal novel mechanisms that are fundamen- 
tal to the ability of cassava to cope with cold stress. 

In summary, our array study will provide the funda- 
mental knowledge related to the biological and physiolo- 
gical changes of cassava under cold stress. It will also be 
served as a very useful genetic resource for relevant 
research community globally. Further studies are needed 
to verify the functions of candidate genes for improving 
cassava tolerance ability through genetic engineering. 

Methods 

Plant materials used for microarray hybridization 

In vitro plants of a cassava cultivar (TMS60444) were 
planted in plastic pots at 28°C under a 16 h light photo- 
period (110-150 (imol-m'^-s'^) for 3 months in the 
greenhouse. Plants with a uniform growth status were 
transferred to a chamber for cold treatment at 7°C 
under weak light (cool-white fluorescent light at 
approximately 35 (imol-m'^-s'^). The apical shoots 
(about 4-6 cm in size. Additional file 1), including apical 
buds, stems, and the first and second expanded leaves, 
were harvested after 4 h and 9 h of treatment, then fro- 
zen in liquid nitrogen and held at -80°C prior to trans- 
porting to Shanghai Biochip Corporation (Shanghai, 
China) for RNA extraction. Untreated apical shoots 
were harvested as controls (0 h). More than 3 plants 
were harvested and pooled for each time point, and the 
collection was repeated 3 times as biological replicates. 

Physiological analyses of cold-treated cassava plants 

To analyze the physiological changes of cassava under 
cold treatment prior to the microarray study, malondial- 
dehyde (MDA) and proline, indicators of the cold 
response in plants, were measured. The MDA content 
was determined by the thiobarbituric acid (TBA) reac- 
tion with minor modifications of the method of Dhindsa 
et al [84] Proline concentrations were measured by the 
sulfosalicylic acid-acid ninhydrin method according to 
Bates et al [85] with slight modifications. 

Determination of hydrogen peroxide (H2O2) content and 
diaminobenzidine (DAB) staining 

Hydrogen peroxide (H2O2) levels were determined 
according to Velikova et al [86] The first and second 
expanded leaves (1 g) were homogenized in an ice bath 
with 10 mL of 0.1% (w: v) trichloroacetic acid (TCA). 



The homogenate was centrifuged at 12,000 x g for 15 
min, and 1 mL of the supernatant was added to 1 mL of 
10 mM potassium phosphate buffer (pH7.0) and 2 mL 
of 1 M KI. The reaction solution containing the super- 
natant was read at 390 nm using a UV-vis spectrophot- 
ometer. The content of H2O2 was given on a standard 
curve. 

H2O2 was visualized by staining with DAB, which 
undergoes a polymerization reaction to yield a dark- 
brown color once it encounters H2O2 [87]. The first or 
second fully developed leaves harvesting from the in 
vitro and greenhouse-grown cassava plants, which were 
treated with cold stress (7°C) for 0, 4, 9, and 24 h, were 
excised. Their leaf petiole cuttings were immersed in 1 
mg/mL of DAB solution (pH3.8) for 8 h under continu- 
ous light. The leaves were immersed in 96% (w/v) boil- 
ing ethanol for 10 min to decolorize the chloroplasts. 
After cooling, the leaves were kept in ethanol and 
photographed. 

Assay for catalase (CAT) and superoxide dismutase (SOD) 
activities 

For the measurement of enzyme activities, the first and 
second fully expanded leaves (1 g) from 3-month-old 
greenhouse plants under cold stress for 0, 4, 9, and 24 h 
were homogenized in 5 mL of 10 mM potassium phos- 
phate buffer (pH7.0) containing 4% (w: v) polyvinylpyr- 
rolidon (Mr 25 000). The homogenate was centrifuged 
at 10,000 rpm for 15 min, and the supernatant obtained 
was used as the enzyme extract. All steps in the pre- 
paration of the enzyme extract were carried out at 0-4° 
C. 

CAT activity was determined by directly measuring 
the decomposition of H2O2 at 240 nm as described by 
Cheng et al [88] The reaction mixture contained 50 
mM potassium phosphate buffer (pH7.0), 10 mM H2O2, 
and 200 (iL of enzyme extract in a 2 mL volume. SOD 
activity assay was based on the method described by 
Beauchamp and Beaucham et al [89], which measures 
the inhibition of the photochemical reduction of nitro 
blue tetrozulium (NBT) at 560 nm. Three milliliters of 
reaction mixture contained 50 mM phosphate buffer 
(pH7.8), 0.1 mM EDTA, 13 mM methionine, 75 (iM 
NBT, 16.7 (iM riboflavin, and 300 (iL of enzyme extract. 

Measurement of total soluble sugars, sucrose and glucose 

Soluble sugars were extracted from the leaf tissues by a 
hot ethanol extraction as follows. Leaves were sampled 
from 3-month-old cassava, which was treated by cold 
stress (7°C) at 0, 4, 9, and 24 h and quickly dried at 
110°C for 15 min and then 70°C overnight. Exactly 100 
mg of homogenized dry leaf powder for each sample 
was extracted with 8 mL of 80% ethanol (v/v) at 85°C 
for 40 min. The extracts were then centrifuged at 
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12,000 X g for 10 min. The ethanol extraction step was 
repeated two times. The 3 resulting supernatants were 
combined and diluted with 80% ethanol to a volume of 
50 mL; they were then treated with 20 mg of activated 
charcoal at 80°C for 30 minutes. The soluble sugar ana- 
lysis was conducted according to Ebell [90] with minor 
modifications. After removing the activated charcoal 
with a 0.2 [im filter, an aliquot of extraction buffer was 
reacted with anthrone-sulfuric acid at 95°C for 15 min 
and then cooled to room temperature. The absorbance 
of the reaction solution was read at 620 nm. The total 
soluble sugar content was calculated by a standard 
curve. Sucrose and glucose content were measured by 
glucose oxidase method as described by Johnson et al 
[91] with slightly modification. 

Oligonucleotide microarray preparation, hybridization, 
and data extraction 

The custom-designed 60-mer microarray [31] was con- 
structed based on public sequence information from a 
large collection of cassava ESTs from NCBI (71,520 
ESTs, released 28 March 2008) and TIGR (Manihot_es- 
culenta_release_5, containing 5,189 assemblies and 
10,214 singletons, released 1 June 2007), as well as a 
35,400 full-length cDNA RIKEN library [27]. Total RNA 
was extracted from both control and cold-treated tissue 
samples using an RNeasy Mini Kit (Qiagen, Valencia, 
CA, USA). The RNA quality was checked on a 1.2% 
agarose gel using an RNase-free electrophoresis system. 
RNA labeling and hybridization were conducted by the 
Shanghai Biochip Corporation (Shanghai, China) follow- 
ing the manufacturer's instructions. The arrays were 
incubated at 65°C for 17 h in Agilent hybridization 
chambers (G2545A) and then washed according to the 
protocol at room temperature. Hybridized microarray 
slides were scanned at 5 \im resolution with an Agilent 
Technologies Scanner (G2505B), and the images were 
saved in JPG format. Both the 10% and 100% Photomul- 
tiplier tube (PMT) settings were selected, and the com- 
bined images were exported. The signal intensities of all 
spots on each image were quantified with the Feature 
Extraction software (Agilent Technologies, Santa Clara, 
California, USA), and the data were saved as .txt files 
for further analysis. 

Data normalization and identification of differentially 
expressed genes (DEGs) 

The signal intensity of each gene was globally normal- 
ized within the GeneSpring GX Software (Agilent Tech- 
nologies, Santa Clara, California, USA) following the 
workflow guide [92]. The signal-to-noise ratio (SNR) 
was calculated as follows: the difference of the median 
signal minus the background median signal, divided by 
the background standard deviation [22]. Pair comparison 



was used to analyze the normalized and averaged data 
from the three types of samples (0, 4, and 9 h). The P 
values were calculated using a t test, and the fold 
changes between each comparison for each gene were 
compared. The genes that were induced or suppressed 
at levels equal to or greater than a twofold ratio (two- 
fold change cutoff, FCC) were taken to be differentially 
expressed when SNR > 2.6 and P value < 0.05. 

Microarray data analysis 

All non-redundant sequences that were considered to be 
unique cassava genes were locally blasted in the TAIR 
protein database (27,217 Arabidopsis protein sequences), 
which was download from ftp://ftp.arabidopsis.org/ 
home/TAIR, using the blastx program in the blastall 
package (version 2.2.9). The top hits were used for gene 
annotations, and the corresponding Arabidopsis gene 
locus identifiers were mapped to the GO function anno- 
tation (http://www.arabiopsis.org/) and to the Kyoto 
Encyclopedia of Genes and Genomes (KEGG) pathways. 
The TAIR percent of cold-responsive genes was calcu- 
lated as follows: The TAIR percent (%) = the number of 
genes annotated to terms in the GO Slim category 
divided by N x 100, where N represents the total num- 
ber of genes from the input list annotated to any term 
in this ontology. 

Utilization of the microarray 

MIAME information about the cassava transcriptome 
microarray used here has been deposited in the Gene 
Expression Omnibus (GEO) of NCBI. The accession 
numbers are: Platform, GPL11271; Series, GSE31073; 
Samples, GSM769563-GSM769571. 

RNA preparation and real-time RT-PCR 

To verify the expression patterns of DEGs identified 
from the microarray analysis, real-time RT-PCR was 
conducted using plant materials from either in vitro 
shoot cultures or greenhouse-grown plants. Four-week- 
old in vitro seedlings were transferred into plastic jars 
containing 150 mL of MS solution [1 x Murashige and 
Skoog basal salts (MS, Duchefa), 1% sucrose, pH5.7] in 
a chamber at 25°C, 70% relative humidity, and a light 
intensity of 125 (imol-m'^- s'^ on a 16 h light/8 h dark 
cycle for 4 days. The solution was then replaced with 
fresh MS medium (pH5.7) supplemented with 100 (iM 
ABA, 25% PEG, or 200 mM NaCl. Meanwhile, the seed- 
lings were directly transferred to a 7°C chamber for cold 
stress treatment. After 6 h, whole seedlings were har- 
vested and frozen in liquid nitrogen; they were then 
stored at -80°C for abiotic stress analysis. Plants in the 
same hydroponic system with no stress treatment were 
harvested and used as controls. The apical buds (AP), 
fibrous roots (FR), mature leaves (ML), stem cambia 
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(SC), young leaves (YL), and young stems (YS) were col- 
lected from 3-month-old cassava plants grown in the 
greenhouse for tissue-specific expression analysis of 
selected genes. All of the stress treatments had three 
biological and temporal replicates. 

The expression of 32 responsive genes was validated 
by real-time RT-PCR with RNA samples extracted with 
a Plant RNA Reagent according to the manufacturer's 
protocol (Invitrogen, Cat.No. 12322-02). The RNA qual- 
ity was determined by running an agarose gel with ethi- 
dium bromide staining. Reverse transcription was 
performed according to the manufacturer's protocol 
(TOYOBO, Code: TRT-101). Each cDNA sample was 
diluted 2 times in sterile ddH20, and 1 (iL of this dilu- 
tion was used as a template for real-time RT-PCR. Spe- 
cific primers were designed with Primer-BLAST (http:// 
www.ncbi.nlm.nih.gov/tools/primer-blast/) to obtain a 
Tm of 60°C and an amplicon length between 70-200 bp 
(Additional file 13). The real-time RT-PCR reactions 
were performed in a 20 \iL volume containing 10 (iL of 
2 X SYBR Green Master Mix (TOYOBO, Code: QPK- 
201), 50 ng cDNA, and 400 nM of forward and reverse 
primers in a Bio-Rad CFX96 thermocycler. The amplifi- 
cation conditions were as follows: 95°C for 1 min, fol- 
lowed by 40-50 cycles of 95°C for 15 s, 60°C for 20 s, 
and 72°C for 20 s. A melting curve (65°C-95°C with a 
heating rate of 0.05°C-s'^ and a continuous fluorescence 
measurement) was run after the PCR cycles. Beta-actin 
was used as the internal control. All of the samples 
were measured in triplicate, and the experiments were 
performed on three biological replicates. The compara- 
tive Ct method was used to calculate the relative gene 
expression level across the samples. The relative expres- 
sion level of each gene in one sample (ACt) was calcu- 
lated as follows: Ct target gene-Ct beta-actin. The 
relative expression of each gene in two different samples 
(AACt) was calculated as follows: ACt (sample 1) -ACt 
(sample 2). 

Additional material 



Additional file 4: Gene descriptions of differentially expressed 
genes based on TAIR protein database. Excel file contains the 319 
DEGs with the highest homologies (AGI locus and full descriptions) and 
their corresponding E values. 

Additional file 5: Gene ontology (GO) annotation of cold responsive 
genes. Excel file contains 266 gene locus identifiers with 370 GO IDs that 
assigned to the three different categories ('Biological Process', 'Molecular 
Function', and 'Cellular Component'). 

Additional file 6: Categories of 'Biological Process' based on GO 
annotation and their contributing genes. 

Additional file 7: Categories of 'Molecular Function' based on GO 
annotation and their contributing genes. 

Additional file 8: Categories of 'Cellular Component' based on GO 
annotation and their contributing genes. 

Additional file 9: GO terms of cold responsive genes related to 
'Response to abiotic and biotic stimulus' and 'Response to stress'. 

Additional file 10: GO terms of cold responsive genes belonging to 
transcription factors. 

Additional file 11: GO terms of cold responsive genes assigned to 
'Chloroplast'. 

Additional file 12: KEGG pathway analysis of cold responsive genes. 

Excel file contains cold responsive genes with mapped KEGG pathways 
and their corresponding pathway entries. 

Additional file 13: Primers used for real-time RT-PCR verification. All 

the forward and reverse primer sequences were included. 
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